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Abstract. The (alpha) version of the Fermi-Pasta-Ulam is revisited through direct 
numerical simulations and an application of weak turbulence theory. The energy 
spectrum, initialized with a large scale excitation, is traced through a series of distinct 
qualitative phases en route to eventual equipartition. Weak turbulence theory is 
applied in an attempt to provide an effective quantitative description of the evolution 
of the energy spectrum. Some scaling predictions are well-confirmed by the numerical 
simulations. 

1. Introduction. One of the very first uses of electronic computing machines was 
Fermi, Pasta, and Ulam's simulation of wave propagation in a weakly nonlinear 



lattice model [16 . They were expecting to observe that the weak coupling of the 
normal modes of the system would induce a redistribution of energy from an initial 
large-scale excitation to an equal sharing (equipartition) of energy among all normal 
modes after some time. As is well known, they were instead surprised to see the 
system display regular behavior characteristic of integrable systems, with the initial 
state recurring on a rather short time scale. This discovery shifted attention to its 
explanation and ramifications Q for several decades. In the last two decades, 
however, the Fermi-Pasta-Ulam (FPU) model has once again been utilized as a test 
model for numerically illustrating and exploring standard concepts in statistical me- 
chanics The peculiar near-integrable behavior observed by Fermi, Pasta, and 
Ulam is characteristic of their model only for systems which are sufficiently small 
in size and energy. There exist well-defined regimes for which the FPU model is 
weakly nonlinear but stochastic p], [24], [35| , and it is in these regimes that one can 
hope to connect the outcome of direct numerical simulations with statistical me- 
chanical concepts JlO|, such as relaxation to equipartition [|l0| pl| |3l], p6[ , entropy 
production |Io|, lij, 15, |3l], chaos as manifested by positive Lyapunov expo- 
nents ||, [iT|lT3|, 35[ 37 1, universal behavior of statistical functions JlJ, |2j| |g[ 3?J 
and virial relations [[3 , f|. The reason for using the FPU model for this purpose is 
that it is one of the simplest and most natural one-dimensional nonlinear models 
for statistical mechanics which can be conceived. An interesting alternative of com- 
parable simplicity which has been the subject of recent research is the truncated 
Burgers model [EL S3, |3(| . 
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Our intent is to use the FPU model to scrutinize a "weak turbulence" (WT) 
theory, a nonequilibrium statistical mechanical theory which attempts to describe 
the dynamical energy transfer among normal modes in a weakly nonlinear, disper- 
sive, extended system [|| [| ^0|, |33|, |34|, The theory has been developed 
over the last four decades to describe the energy transfer in wave dynamics pri- 
marily in fluids and plasmas [33| |34| , f4(i| ], among other novel applications such as 
semiconductors |2^, |2(], |27]] . 

These systems are generally too complex for effective comparison of the weak 
turbulence theory with direct numerical simulations. Only in recent years has a 
simple one-dimensional model with features representative of such fluid systems 
been explored by Cai, Majda, McLaughlin, and Tabak |Q |[ to examine the 
assumptions underlying WT theory. We propose to apply the FPU model for a 
similar purpose, though the issues on which we focus are distinct. The implications 
of our studies for the framework of WT theory will be taken up in other works pl[| . 

Here we will present the content of our findings as they inform the relaxation 
process in the FPU model in the stochastic but weakly nonlinear regime. We 
will restrict attention to the a version of FPU model, which has purely quadratic 
nonlinearity in the equations of motion (Section [2]). The (3 version (with cubic 
nonlinearity) seems to be the subject of more work j|, |l5|, [22| |3f], [55| |36|], but 
the a version has attributes which make it more suitable for a first test case for 
WT theory. Most of the previous statistical mechanical work concerning the FPU 
models of which we are aware focuses primarily on computing particular statistical 
measures of the process, such as the time until equipartition is reached [3q] , 
the Lyapunov exponents characterizing the degree of chaos 0, |13|, pi p| , or more 
exotic quantities characterizing the geometry of the trajectories lf|. Another 
recent line of research has been tracing the path of energy transfer starting from a 
small set of excited modes [58], p9| . Because the WT theory has the potential 
to describe the process of energy transfer in the system from beginning to end, 
we have instead sought to characterize the entire evolution of the energy spectrum 
from large-scale excitation to eventual equipartition. We will consider the energy 
transfer in spectral terms, in contrast to the physical space viewpoint developed for 
the /3-FPU model by Lichtenberg and coworkers fjlf . 

The energy spectrum in the a-FPU model approaches equilibrium through a se- 
ries of qualitatively distinct phases which we illustrate in Section At the initial 
time, the energy is concentrated in a small set of low-wavenumber modes. This en- 
ergy then proceeds to higher wavenumbers first through a standard superharmonic 
cascade, and then shifts to a nonlocal transfer of energy from low wavenumbers 
to a band of intermediate wavenumber modes. The energy in this intermediate 
wavenumber band then rolls back through an inverse cascade to lower wavenum- 
bers again. This process then creates approximate equipartition only over a set of 
modes extending up to a cutoff wavenumber, beyond which the energy content falls 
off exponentially rapidly |9|, |T^] . We refer to the location of the transition between 
the flat and rapidly decaying parts of the energy spectrum as the "knee." 

After presenting this pictorial "life history" of a large-scale excitation in the a- 
FPU system, we present in Section [| some specific quantitative predictions of WT 
theory and compare them with the numerical results. At the coarsest level, WT 
theory suggests the presence of two nonlinear time scales. Over the first nonlinear 
time scale, energy is exchanged through triads of modes which remain resonant over 
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this time scale. A consideration of the resonances in the dispersion relation indi- 
cates that only modes of sufficiently small wavenumber can participate in nearly 
resonant triads ]37| . This in turn suggests that this triad interaction phase should 
correspond to the formation of the partial equipartition up to the knee. The subse- 
quent relaxation of the energy spectrum to global equipartition requires the slower 
energy exchange among resonant quartets of normal modes, which are much more 
abundant. Our present focus is on the triad interaction phase. 

An adaptation of the WT theory allows predictions of both the order of mag- 
nitude of the time scale and the location of the knee which agree excellently with 
direct numerical simulation. 



2. Fermi-Pasta-Ulam Model. The Fermi-Pasta-Ulam (FPU) model is a model 
for a one-dimensional collection of particles with massless, weakly anharmonic (non- 
linear) springs connecting them to each other. Letting {qj} jLi and {pj}f =1 denote 
the position and momentum coordinates of an A-particle chain, we can define the 
FPU model Hamiltonian: 

H = y ( *L + >--tt+i) 2 + a fa-tt+i) 3 + ^-i^r) (1) 

r— J y2m 2 3 4 J 

Here we assumed periodic boundary conditions pn+i = Pi and qN+i — qi- Equiva- 
lently, the beads are connected in a circular arrangement. The parameter m denotes 
the particle mass, while k, a, and (3 are coefficients involving the spring properties. 
The equations of motion are the standard Hamilton's equations: 

• dH . 3H 

* = WS Pj = ~Wj (2) 

We will here focus on the a-FPU model for which a / and the quartic term is 
absent ((3 = 0). We nondimcnsionalize the system with respect to the spring con- 
stant k, the mass m, and the energy density H/N. Retaining the original symbols 
for the nondimensionalized variables pj and qj, we obtain the nondimcnsionalized 
Hamiltonian and equations of motion: 

N N 

H = 2£(P? + fo-W) + f£fo-fc-+i) 8 ( 3 ) 
Qj = Pj 

Pj = (Qj-i - 2 * +q 3 +i) (1 +efe-i -Qj+i)) 
Our choice of nondimensionalization implies that 

H = N (4) 

for all times. The fundamental nondimensional parameter measuring the strength 
of the nonlinearity is 
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In order to study the transfer of energy among different scales, we represent the 
system in terms of Fourier modes: 
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The Hamiltonian in the new variables reads 

N_ 

n = ^ E m 2 +^i\Q k \ 2 ] 



(5) 



(6) 



fe=-f +i 



37V 2 



22 M ,ka Qki Qk 2 Qk 3 S kl +k 2 + k 3 ,0 , 



where the dispersion relation is given by 



LU k = 2 



sm 



irk 



the nonlinear coupling coefficients are 

V ku k 2 ,k 3 = -isgn(kik 2 k 3 )uj kl uj k2 uJk a , h + k 2 + k 3 = 0, 

and 

r j 1 if i = j mod N, 
else. 



(7) 
(8) 



J i,3 



is a periodized version of the Kronecker delta function. 

To quantify the amplitude of activity of the FPU chain at different scales, we 
define the harmonic energy contribution of each Fourier mode: 

Eh,k(t) = ^[\Pk\ 2 +^k\Qk\ 2 } ■ 

Energy equipartition implies Eh, k is independent of k (and t). The total harmonic 
contribution to the energy is 

N/2 

E h (t) = E 'M- 

k=- N/2+l 

3. Numerical Simulation of Relaxation from Initial Large-Scale Excita- 
tion. The equations of motion in (^) were integrated using a fourth order Runge- 
Kutta time stepping routine. Since the equations are local in physical space, all 
calculations were performed there. FFTs were employed only for data analysis. All 
simulations conserved energy to within 10~ 3 after the full length of the simulation 
whereas total momentum and particle position were conserved to machine accuracy. 
With a specified number of excited initial modes, the amplitude and phase of the 
Fourier coefficients are each chosen from a uniform distribution so that, on aver- 
age, each mode would be initialized with the same energy. The initial data is then 
normalized according to Eq. M). The results of all experiments are averaged over 
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10-20 independent realizations of the initial data; these relatively small ensembles 
proved sufficient to elucidate the results. 

Figure [l] shows a sequence of ensemble averaged spectra for a lattice of length 
N = 512 and nonlinearity strength e = 0.03 at times t = 0, 50, 100, 200, 400, 
1000, 2500, 5000, 10 4 , displaced for ease of viewing. The evolution proceeds from 
a set of 50 initially excited modes to a superharmonic cascade to all wave numbers 
with exponentially decreasing energy by t = 50. By t = 100 the initial band has 
transferred much of its energy to intermediate wavenumbers, forming a slight hump. 
Thereafter, this hump of energy rolls back via an inverse cascade to low wavenum- 
bers. At t = 10 4 , the energy spectrum exhibits a plateau at low wavenumbers and 
an exponential falloff at higher wavenumbers. This last spectrum is the motivation 
for the term "knee", below which the waves are in equipartition and above which 
they are not substantially excited ||, [l?]]. After t — 10 4 the spectrum evolves over 
much longer time scales, eventually arriving at equipartition throughout. 

Figure ^ shows a similar experiment with a larger ensemble (20 realizations), 
N = 512 and e = 0.05. The initial excitation band is very much smaller, including 
only 20 waves. Again spectra are displaced and in this example are shown at 
t = 0, 10 3 ,2 x 10 3 ,4 x 10 3 ,8 x 10 3 ,2 x 10 4 ,5 x 10 4 , 10 5 . Energy is driven first 
to an intermediate range of wavenumbers which saturate. Subsequently, an inverse 
cascade of energy extends the band of equilibration backward to lower wavenumbers 
until t ~ 8 x 10 3 , at which point only the lowest wavenumber has yet to reach 
equipartition. At t = 5 x 10 4 the spectrum is quasi-stationary, equipartition being 
achieved among all wavenumbers less than fcknco- The energy in larger wavenumbers 
decreases rapidly. At t = 10 5 the highest modes begin to acquire energy. Eventually 
the whole spectrum will arrive at equipartition; this process is outside of the scope 
of the current work. 

4. Scaling Predictions. A renormalized WT theory derived in [^lj predicts that 
significant three wave interaction should occur in a band 

|fc| < fckncc ~ N^~e (9) 

on a time scale 

T 3 ~ e~ 3/2 . (10) 
This theory only applies when the lattice size is large enough (N e -1 / 2 ) and the 
number of initially excited modes is an order unity fraction of the knee width fcknee, 
so that the renormalized energy spectrum remains self-consistently of order unity 
during this phase of evolution. 

A useful statistical measure for our purposes is the spectral entropy, defined as 

This provides a measure of the effective number of excited normal modes at any 
given time, n cS (t) = e s M gfj g g |, Figure | shows rescaled plots of 
this spectral entropy as a function of time. The onset of the quasi-stationary 
phase, after the end of the three-wave evolution, is clearly evident. The knee width 
fcknee ~ l-5Ny/e is determined as an average of n e ff over a time window shortly after 
the entropy ends its rapid rise. This scaling relationship is robust against various 
choices of initial bandwidth excitations. The time to reach partial equipartition T3, 
however, does depend sensitively on the choice of initial data. As discussed above, 
the WT theory producing the scaling prediction ( [To| ) assumes the initial data is 
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excited over a band of wavenumbers which is an order unity fraction of the knee 
width. To test the prediction (|i"o|), then, we choose the system to have initially 
ifcknee ~ 0.75N^/e excited modes. (The evolution depicted in Figure |l| comes from 
initial data of this form, whereas the example in Figure || was initialized with a 
considerably smaller band of excited modes). The time scale T3 for the system to 
reach partial equipartition is determined automatically as the first time at which 
n e s{t) achieves the value &kncc = 1.5Ny/e. 

4.1. Effect of strength of nonlinearity. Ensembles of FPU lattices with lattice 
length N — 512 were integrated for values of e ranging from 10~ 3 to lO^ 1 . Figures 
[| and H show excellent agreement with the e scaling dependences (|]) and ( |io| ) 
predicted by the renormalized WT theory. The scaling fck nC c ~ e 1 ^ 2 had also been 
previously observed in |^7| . 

4.2. Effect of lattice size. The renormalized WT theory also predicts that the 
knee width should scale with lattice size (Eq. ([)])) whereas the three wave time 
scale should not (Eq. (|l0|)). These properties are compatible with a thermodynamic 
limit. In figure || the logarithm of the fraction of modes to the total lattice size is 
plotted against time for seven ensembles of experiments with increasing lattice size 
N = 32,64,128,256,1024,2048,4096, and fixed e = 0.1. These experiments were 
again initialized with half the number of modes of the predicted knee. Three features 
of this plot stand out most clearly. First, the spectral entropy follows a universal 
evolution |lj, ||(| for lattice sizes larger than N = 128. Secondly, the number of 
excited modes exponentially increases with time prior to three wave equilibrium. 
Finally, for N — 32 (where the initial number of excited modes is approximately 
7) there is no equilibrium, but rather quasi-periodic behavior. In fact, a shadow 
of this behavior is present for TV = 64 and 128 also. These are reminiscent of the 
integrability discovered in Fermi, Pasta and Ulam's original work jl6| , and indicate 
the breakdown of the WT theory scaling predictions for small lattice sizes. 

5. Conclusion. We have emphasized some of the long-lived transient features of 
the evolution of the energy spectrum in the a-FPU model. Weak turbulence theory 
has been successful in predicting scaling exponents concerning the achievement of 
partial equipartition. In the future, we will endeavor to explain other dynamical 
aspects, such as the formation of the energy hump at intermediate wavenumbers 
and subsequent inverse cascade, in similar quantitative terms. 

6. Acknowledgements. The authors would like to thank David Cai and Gregor 
Kovacic for helpful discussions. JAB is supported by an NSF VIGRE postdoctoral 
research fellowship DMS 9983646, PRK is supported by an NSF grant DMS-A11271, 
and YVL is partially supported by NSF Career grant DMS 0134955 and ONR YIP 
grant N000140210528. 

REFERENCES 

[1] Rafail V. Abramov, Gregor Kovacic, and Andrew J. Majda. Hamiltonian structure and sta- 
tistically relevant conserved quantities, to appear in Communications on Pure and Applied 
Mathematics, 2002. 

[2] Carlo Alabiso and Mario Casartclli. Quasi-harmonicity and power spectra in the FPU model. 

J. Phys. A, 33(5):831-839, 2000. 
[3] Carlo Alabiso and Mario Casartelli. Normal modes on average for purely stochastic systems. 

J. Phys. A, 34:1223-1230, 2001. 



STAGES OF ENERGY TRANSFER IN THE FPU MODEL 



7 




lo 9ioK! 



Figure 1. Temporal evolution of spatial FPU energy spectrum 
versus mode frequency for an ensemble of experiments on a lat- 
tice of length N = 512. Subsequent spectra are shifted up- 
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32,64,128,256,1024,2048,4096, e = 0.1. The simulations with 
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show some mild quasi-periodicity whereas the simulation N = 32 
clearly shows the quasi-periodic behavior of the original FPU sim- 
ulations ||l6| . 
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